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ABSTRACT 


A  Monte  Carlo  technique  is  described  for  simulating  the 
high-speed  flow  of  a  rarefied  gas  around  a  cylindrical  body. 

The  method  uses  a  hard-sphere  model,  or,  alternatively,  the 
"Maxwellian"  model  to  compute  collisions  between  particles.  A 
representative  "test  particle"  is  scattered  in  such  a  way  that 
its  statistical  behavior  approximates  that  of  a  random  particle 
chosen  from  the  actual  gas.  This  requires  knowledge  of  the 
distribution  function  of  "target"  particles,  which  is  estimated 
from  observations  of  the  behavior  of  the  test  particle. 

The  method  is  applicable  to  any  collision  model  and  to  any 
model  for  the  interaction  between  the  gas  and  the  surface  of 
the  body. 


Accepted  for  the  Air  Force 

Franklin  C.  Hudson 

Chief,  Lincoln  Laboratory  Office 


iii 


TABLE  OF  CONTENTS 


ABSTRACT  Hi 

I.  INTRODUCTION  1 

II.  REVIEW  OF  PREVIOUS  MONTE  CARLO  GAS  SIMULATIONS  3 

III.  BRIEF  DESCRIPTION  OF  THE  PRESENT  CALCULATION  8 

IV .  MATHEMATICAL  BASIS  OF  THE  SIMULATION  18 

(1)  Dynamics  of  elastic  collisions  18 

(2)  Collision  mechanics  for  a  general  repulsive  19 

force  law 

(3)  Maxwellian  molecule  22 

(4)  Hard  sphere  molecule  24 

(5)  Calculation  of  velocities  after  collision  24 

(6)  Selection  of  collision  parameters  25 

(7)  Calculation  of  collision  rates  27 

(8)  Mean  free  path  30 

(9)  Representation  of  the  distribution  function  31 

of  a  finite  collection  of  particles 

(10)  Deciding  when  the  next  collision  will  occur  34 

(11)  Choice  of  a  collision  partner  35 

(12)  Estimation  of  the  distribution  function  36 

(13)  Deletion  of  entries  from  the  lists  of  delta  40 

functions 

(14)  Boundary  conditions  40 

(15)  Generation  of  new  test  particles  41 


v 


(16)  Interaction  of  particles  with  the  aerodynamic  45 
surface 

(17)  Sampling  48 
V.  CONCLUSIONS  51 
REFERENCES  52 
APPENDIX  53 


vi 


I.  INTRODUCTION 

Many  attempts  have  been  made  to  solve  analytically  the 
problem  of  high  speed  flow  of  a  rarefied  gas  around  various 
simple  aerodynamic  bodies.  Although  some  progress  has  been 
made,  exact  solutions  are  not  yet  available  even  for  very  sim¬ 
plified  boundary  conditions  and  gas  models.  When  the  physics 
of  the  problem  is  described  more  realistically,  it  becomes  very 
difficult  to  get  even  approximate  solutions  which  agree  with 
experimental  observations. 

One  technique  which  may  be  used  to  circumvent  the  analyt¬ 
ical  difficulties  is  the  "direct"  Monte  Carlo  simulation.  In 
such  an  approach,  the  gas  is  modeled  by  a  collection  of 
"particles",  generally  much  fewer  in  number  than  the  number  of 
particles  in  the  actual  gas.  These  particles  are  made  to  be¬ 
have  as  actual  gas  particles  would,  that  is,  they  are  affected 
by  the  same  statistical  laws  that  govern  the  actual  distribution 
of  particle  velocities  and  positions.  By  observing  these  "test" 
particles  for  a  sufficiently  long  time,  one  can  then  deduce 
their  distribution  function,  which  when  properly  normalized 
should  be  the  distribution  function  for  particles  in  the  actual 
gas.  Since  all  aerodynamic  properties  can  be  deduced  from  the 
distribution  function,  it  constitutes  a  complete  solution  to 
the  problem. 

An  "incident"  particle  can  be  made  to  behave  statistically 
as  if  it  were  a  typical  member  of  the  actual  gas  only  if  one 
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knows : 

(1)  the  mechanism  governing  collisions  between  particles, 

(2)  the  interaction  between  gas  particles  and  any  solid 
surface  which  they  may  strike, 

(3)  the  density  and  velocity  distribution  of  the  "target" 
particles  which  scatter  the  incident  particle. 

The  first  requirement  is  easily  met  by  postulating  a 
collision  model,  either  based  upon  empirical  data  or  upon  some 
repulsive  force  law  chosen  for  analytic  simplicity.  Similarly, 
a  model  for  the  particle-surface  interaction  may  be  chosen. 

The  third  requirement,  however,  requires  knowledge  of  the  dis¬ 
tribution  function,  which  is  the  object  of  the  calculation. 

Thus,  the  simulation  envisioned  here  must  be  iterative  in  nature, 
hopefully  converging  to  the  correct  answer  after  a  sufficiently 
long  run  under  constant  conditions.  The  test  particles  must  be 
scattered  by  a  collection  of  particles  distributed  in  velocity 
space  according  to  the  latest  estimate  of  the  true  distribution. 
In  turn,  it  must  be  possible  to  improve  this  estimate  by  incor¬ 
porating  information  gained  by  observing  the  behavior  of  the 
scattered  particles. 
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II.  REVIEW  OF  PREVIOUS  MONTE  CARLO  GAS  SIMULATIONS 

Only  a  few  attempts  have  been  made  at  solving  problems  in 

gas  dynamics  by  direct  simulation.  One  of  the  earliest  was  by 

Alder  and  Wainwright'*" ,  who  followed  the  exact  motion  of  about 

2 

100  rigid  elastic  particles.  Bird  solved  by  a  Monte  Carlo 
method  for  the  relaxation  to  thermal  equilibrium  of  about  1000 
hard  sphere  particles,  starting  with  all  the  particles  moving 
at  the  same  speed  but  in  random  directions.  A  similar  calcu¬ 
lation,  using  two  different  collision  models,  has  recently  been 
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performed  here. 

The  calculations  described  above  were  essentially  zero¬ 
dimensional  (in  position  space),  since  all  the  particles  were 

in  effect  located  at  a  single  point. 

4  5  6 

Lavin  and  Haviland  ’  ’  treated  two  more  significant 

problems  by  the  Monte  Carlo  method:  heat  transfer  between  two 
parallel  planes,  and  the  one-dimensional  shock  wave.  Both 
problems  are  one-dimensional  in  position  space,  and  symmetry 
allowed  the  problem  to  be  posed  in  a  three-dimensional  phase 
space  consisting  of  the  position  along  the  line  normal  to  the 
planes,  the  velocity  along  this  same  line,  and  the  magnitude 
of  the  velocity  component  parallel  to  the  planes.  With  only 
three  dimensions,  it  was  feasible  to  divide  all  of  phase  space 
into  cells  of  appropriate  size.  In  each  cell  was  stored  a 
single  scalar  number  giving  the  value  of  the  distribution  func¬ 
tion.  A  single  test  particle  was  introduced  at  one  surface 


3 


with  a  random  thermal  velocity  corresponding  to  the  temperature 
at  that  surface.  It  was  then  advanced  ballistically  from  cell 
to  cell  in  phase  space.  If  it  struck  a  surface,  it  was  re¬ 
emitted  randomly  at  the  temperature  of  the  surface.  When  neces¬ 
sary,  collisions  were  computed.  This  could  easily  be  done, 
since  an  approximation  to  the  distribution  function  was  always 
available.  By  observing  the  single  test  particle  for  a  long 
time,  the  current  approximation  to  the  distribution  function 
could  then  be  improved.  This  method  attempts  an  exact  solution 
to  the  problem  (i.e.,  the  distribution  at  all  points  in  space) 
and  is  limited  by  the  coarseness  of  the  grid  of  cells  and  by 
the  computing  time  required  for  statistical  fluctuations  to 
average  out.  It  is  also  evident  that  if  the  Knudsen  number  (the 
ratio  of  mean  free  path  to  some  characteristic  length)  were  too 
small,  an  excessive  number  of  collisions  would  have  to  be  com¬ 
puted.  Finally,  one  would  prefer  an  initial  distribution  which 
is  somewhere  near  the  expected  true  distribution.  Otherwise,  it 
is  not  clear  that  the  test  particle  would  be  scattered  in  such 
a  way  as  to  yield  an  improved  estimate. 

While  the  method  of  Haviland  and  Lavin  produces  good  results 
for  one-dimensional  problems,  it  is  not  easily  extended  to  prob¬ 
lems  lacking  the  symmetry  needed  to  reduce  phase  space  to  three 
dimensions.  To  divide  a  five-  or  six-dimensional  phase  space 
into  a  grid  of  cells,  even  if  the  grid  is  very  coarse,  would 
exceed  the  storage  capacity  of  most  computers.  Furthermore, 
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the  distribution  function  (which  is  proportional  to  the  frac¬ 
tion  of  the  time  a  given  particle  spends  in  a  particular  cell 
in  phase  space)  is  estimated  by  averaging  the  behavior  of  the 
test  particle  for  a  very  long  time.  Thus  the  test  particle 
must  enter  each  cell  a  large  number  of  times  before  meaningful 
averages  can  be  obtained.  Thus  not  only  the  storage  require¬ 
ments  but  also  the  computation  time  increases  by  many  orders 
of  magnitude  in  going  from  three-  to  five-  or  six-dimensional 
phase  space. 

Fortunately,  in  most  gas  dynamics  problems,  only  certain 
averages  of  the  distribution  function  over  all  velocity  space 
are  required,  e.g.,  density,  temperature,  mean  velocity,  etc., 
rather  than  the  velocity  distribution  itself.  Such  averages  can 
be  obtained  quite  accurately  even  from  a  rather  crude  approxi¬ 
mation  to  the  distribution  function,  provided  that  the  approxi¬ 
mation  is  not  biased.  Thus  a  Monte  Carlo  method  may  provide 
useful  information  even  when  there  is  not  enough  symmetry  to 
reduce  phase  space  to  three  dimensions. 

Probably  the  greatest  success  in  this  direction  has  been 
7 

achieved  by  Bird  ,  who  applied  a  novel  method  to  problems  with 
either  axial  or  cylindrical  symmetry.  In  these  problems  (which 
include  gas  flow  around  cylinders,  spheres,  flat  plates,  cones, 
etc.)  phase  space  is  reduced  to  five  dimensions  by  ignoring  one 
position  coordinate.  Bird  stores  the  two  positions  and  three 
velocity  components  of  a  small  number  of  particles  (about  1000) 


5 


to  represent  his  model  gas,  and  does  not  divide  phase  space  into 
cells.  In  order  to  determine  whether  or  not  a  given  "incident" 
particle  collides  at  a  particular  location,  one  must  have  some 
representation  of  the  distribution  function  of  "target"  particles 
at  that  point.  Bird  uses  as  this  representation  the  particles 
which,  at  the  given  instant,  happen  to  lie  within  a  sphere  cen¬ 
tered  at  the  would-be  collision  site.  This  sphere  is  large 
enough  to  give  a  reasonably  large  choice  of  collision  partners. 

A  large  number  of  such  collisions  are  calculated  (the  details 
of  this  process  are  described  in  Bird's  paper)  using  the  "hard 
sphere"  collision  model  with  the  entire  collection  of  particles 
frozen  in  position.  Then  all  the  particles  are  advanced  ballis- 
tically  (without  collisions)  for  a  time  equal  to  the  accumulated 
mean  waiting  times  for  each  of  the  collisions  which  were  calcu¬ 
lated.  The  procedure  is  done  in  such  a  way  that  each  particle 
has  the  prescribed  mean  free  path.  Samplings  are  made  at  inter¬ 
vals  long  enough  to  permit  decorrelation,  and  the  desired  aero¬ 
dynamic  quantities  are  obtained  from  ensemble  averages  of  many 
samples.  Bird  presents  density  profiles,  drag  coefficients, 
and  various  other  aerodynamic  quantities  of  interest  for  various 
Knudsen  numbers,  stream  speeds,  body  temperatures,  and  body 
geometries . 

While  Bird's  method  produces  impressive  results,  the  proc¬ 
ess  by  which  convergence  occurs  is  a  subtle  one  whose  validity 
is  not  easy  to  eastablish.  Thus  it  seems  desirable  to  find  some 
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alternative  Monte  Carlo  technique  whose  results  could  be  com¬ 
pared  with  Bird's.  Such  a  method  will  be  described  in  detail 
below. 
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III.  BRIEF  DESCRIPTION  OF  THE  PRESENT  CALCULATION 


The  model  used  in  this  calculation  may  be  regarded  as  a 
compromise  between  the  technique  of  Haviland  and  Lavin  and 
Bird’s  method.  Haviland  and  Lavin  divided  all  of  phase  space 
into  cells.  Bird  used  no  cells  but  rather  a  single  list  of 
particles  to  represent  his  distribution  function.  For  this 
calculation,  position  space  is  divided  into  a  two-dimensional 
grid  of  cells,  while  the  velocity  distribution  in  each  cell 
is  represented  by  a  list  of  four-vectors.  The  first  three  com¬ 
ponents  are  the  velocity  components  of  a  "target  particle", 
while  the  fourth  component  is  a  weighting  factor. 

In  a  given  cell,  the  true  continuous  velocity  distribution 
is  thus  replaced  by  a  collection  of  delta  functions  whose  loca¬ 
tions  are  given  by  the  first  three  components  and  whose  areas 
are  given  by  the  fourth  component.  When  calculations  of  colli¬ 
sion  rates  and  selection  of  collision  partners  call  for  a  veloc¬ 
ity  distribution,  this  collection  of  "spikes"  is  normalized  by 
dividing  by  the  sum  of  all  the  fourth  components.  It  may  then 
be  treated  mathematically  as  if  it  were  the  true  continuous 
velocity  distribution.  The "particles"  on  this  list  (a  different 
list  for  each  cell)  have  no  purpose  other  than  to  serve  as  a 

local  representation  of  the  velocity  distribution.  They  are 
not  moved  ballistically  from  cell  to  cell,  as  are  Bird’s  par¬ 
ticles  . 
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Figure  1.  The  grid  of  cells  in  position  space.  The  radius  of  the 
cylinder  is  /45.  The  lower  left  corners  of  some  typical 
cells  are  listed  below,  in  order  to  illustrate  the  num¬ 
bering  system. 
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0) 
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24) 
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at 
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(67,  12) 

at 

(  24, 
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There  are  a  total  of  804  cells. 
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The  formation  of  the  lists  and  their  statistical  justi¬ 
fication  will  be  described  later.  The  actual  grid  of  cells 
used  is  shown  in  Figure  1.  Because  finer  structure  and  higher 
density  are  expected  near  the  cylinder  than  far  away  from  it, 
three  different  cell  sizes  are  used  in  the  grid.  The  stream 
velocity  is  in  the  x-direction*  and  the  axis  of  the  cylinder  is 
in  the  z-direction.  The  entire  grid  is  bounded  by  the  plane  of 
symmetry  (y  =  0)  and  by  upstream,  downstream,  and  side  boundaries 
far  from  the  cylinder.  It  is  assumed  that  collisions  outside 
these  boundaries  will  not  appreciably  affect  conditions  near 
the  cylinder. 

Like  the  calculation  of  Haviland  and  Lavin,  this  method 
used  only  one  "test  particle"  which  wanders  from  cell  to  cell. 

The  process  begins  by  generating  a  test  particle  at  a  random 
position  along  an  end  or  side  boundary  with  the  free  upstream 
velocity  plus  a  randomly  generated  thermal  velocity. 

The  test  particle  then  enters  the  first  cell,  and  the  time 
required  to  traverse  the  cell  is  computed.  The  collision  rate 
(which  depends  on  the  collision  model,  the  test  particle  veloc¬ 
ity,  the  local  density  and  the  local  velocity  distribution)  is 
calculated,  and  a  random  "time-to-collision"  is  then  drawn  from 


A  component  of  stream  velocity  in  the  z-direction  could  also 
be  added  without  destroying  symmetry,  but  this  has  not  been 
done  in  the  first  calculations. 
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the  appropriate  distribution.  If  this  time  is  greater  than 
the  time  to  traverse  the  cell,  no  collision  occurs  and  the  test 
particle  is  advanced  to  the  beginning  of  the  next  cell.  (The 
traversal  of  the  cell  by  the  test  particle  clearly  provides  in¬ 
formation  about  the  distribution  function;  this  information  is 
incorporated  into  the  current  estimate  for  the  cell  by  a  process 
to  be  described  later.)  If  the  time  to  collision  is  less  than 
the  time  to  the  cell  wall,  a  collision  will  occur.  The  test 
particle  is  then  first  advanced  ballistically  to  the  site  of  the 
collision,  and  the  information  gained  from  this  portion  of  the 
trajectory  is  incorporated  into  the  current  estimate  for  the 
cell.  Then  a  velocity  for  the  "collision  partner"  is  drawn 
randomly  from  the  list  (the  actual  drawing  depends  on  the  colli¬ 
sion  model)  and  a  pair  of  collision  parameters  (again  dependent 
on  the  collision  model)  is  drawn.  The  new  velocity  of  the  test 
particle  is  then  calculated.  The  list  of  velocities  representing 
the  distribution  function  is  not  directly  affected  by  the  colli¬ 
sion. 

Using  the  new  velocity,  the  time  to  the  wall  and  the  time 
to  the  next  collision  are  recalculated,  and  the  process  continues. 
There  is  no  inherent  limitation  on  the  number  of  collisions  which 
may  occur  during  one  passage  through  a  cell,  but  multiple  colli¬ 
sions  during  one  passage  will  be  rare  due  to  the  long  mean  free 
path.  The  number  of  updatings  (the  incorporation  of  information 
from  a  straight  segment  of  the  test  particle  trajectory) 
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corresponding  to  a  single  passage  through  a  cell  is  one  greater 
than  the  number  of  collisions  occurring  during  that  passage. 

Thus  the  test  particle  wanders  from  cell  to  cell,  occasion¬ 
ally  undergoing  a  collision  which  changes  its  velocity.  Whether 
or  not  it  collides  in  a  particular  cell,  the  fact  that  it  passed 
through  contains  information  about  the  local  density  and  veloc¬ 
ity  distribution.  This  information  must  be  added  to  the  current 
local  estimate  to  produce  an  improved  estimate.  One  way  to  do 
this  would  be  simply  to  add  one  more  velocity  to  the  list, 
namely  the  test  particle  velocity.  The  weighting  factor  would 
be  the  time  the  particle  spent  in  the  cell  during  its  passage. 
This  is  necessary  in  order  that  particles  which  remain  in  the 
cell  a  long  time  contribute  more  to  the  density  than  those  par¬ 
ticles  which  go  through  very  quickly.  This  also  takes  proper 
account  of  trajectories  which  barely  cut  across  the  corner  of 
the  cell. 

This  procedure  would  eventually  produce  lists  which  would 
be  too  long  for  the  storage  capacity  of  the  computer.  Even  if 
this  were  no  problem,  the  distribution  function  represented  by 
a  very  long  list  would  not  be  appreciably  changed  by  the  addi¬ 
tion  of  one  more  four-vector,  so  convergence  would  be  very  slow. 
Thus,  the  actual  procedure  followed  is  to  let  the  list  grow  to 
a  maximum  length  (not  more  than  22  in  this  version).  After  the 
maximum  length  has  been  reached,  an  entry  is  deleted  every  time 
a  new  entry  is  added.  Thus  an  entry  remains  on  the  list  for 
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n  passages  of  the  test  particle,  where  n  is  the  length  of  the 
list . * 

While  the  deletion  of  entries  from  the  list  destroys  some 
of  the  information,  the  features  of  the  distribution  function 
of  greatest  interest,  e.g.,  density,  mean  velocity,  temperature 
and  certain  other  higher  moments  of  the  distribution,  are  accu¬ 
mulated  in  positions  on  the  list  reserved  for  this  purpose. 

Thus  the  current  estimates  of  mean  velocity  and  density  are  not 
based  merely  upon  the  small  collection  of  four-vectors  currently 
on  the  list,  but  upon  all  entries  which  have  been  on  the  list 
since  the  last  time  the  velocity  and  density  accumulators  were 
reset  to  zero. 

When  the  test  particle  strikes  the  body,  it  is  absorbed 
and  re-emitted  at  the  same  point  according  to  the  model  chosen 
to  represent  the  surface.  A  procedure  often  used  for  analytic 
work,  and  also  used  here,  is  to  re-emit  a  fraction  a  diffusely, 
and  the  rest  specularly.  Diffuse  emission  is  equivalent  to 
having  a  small  hole  in  the  surface  at  the  absorption  point, 
behind  which  is  a  gas  in  equilibrium  at  some  temperature  T^. 

When  a  particle  is  absorbed,  a  shutter  in  front  of  the  hole  is 
opened  and  a  single  particle  is  allowed  to  pass  outward.  This 
sort  of  emission  follows  the  usual  cosine  law,  and  the  velocity 


5|t 

An  alternative  representation  of  the  velocity  distribution  is 
described  in  the  appendix. 
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of  the  emitted  particle  is  uncorrelated  with  that  of  the  ab¬ 
sorbed  particle.  Specular  reflection  simply  reverses  the  normal 
component  of  velocity.  The  fraction  a  which  is  emitted  diffusely 
is  called  the  accommodation  coefficient. 

This  surface  model,  with  a  near  unity,  is  considered  rea¬ 
sonable  for  the  low  densities  of  interest  in  this  problem.  How¬ 
ever,  it  would  not  be  difficult  to  use  a  more  complicated  model 
based  on  experimental  observations. 

When  the  test  particle  strikes  the  plane  of  symmetry,  it 
is  reflected  specularly.  When  it  crosses  a  side  or  end  boundary, 
it  is  lost.  A  new  test  particle  is  then  generated  with  position 
and  velocity  uncorrelated  with  those  of  the  old  test  particle. 

A  certain  fraction  of  the  particles  must  be  generated  along  the 
side  boundary  in  order  to  compensate  for  downstream  depletion 
through  the  sides  (which  would  occur  even  if  the  body  were  ab¬ 
sent)  .  This  fraction  depends  upon  the  ratio  of  the  stream  speed 
to  the  thermal  speed  and  is  calculated  at  the  beginning  of  the 
computation. 

The  velocity  of  the  newly  generated  test  particle  is  then 
formed  by  adding  a  random  thermal  velocity  to  the  x-directed 
stream  velocity.  If  a  particle  is  to  be  generated  at  an  end 
boundary  (a  result  of  drawing  a  rectangular  random  number  greater 
than  the  fraction  to  be  introduced  at  the  side),  it  is  given  a 
random  y-coordinate  along  the  upstream  or  downstream  boundary, 
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depending  on  the  sign  of  the  x-component  of  velocity.  If  it 
is  introduced  at  the  side  boundary,  its  y-velocity  is  always 
negative . 

After  a  predetermined  number  of  particles  have  been  lost 
through  the  end  and  side  boundaries,  a  sample  is  taken.  Only 
the  "cumulative  registers"  in  the  individual  cells  contribute 
to  the  sampling.  From  these  cumulative  registers  are  formed 
arrays  showing  the  density  (normalized  to  free-stream  density), 
the  mean  velocity,  the  thermal  velocity  spread  and  the  energy 
flow  vector  (normalized  to  free-stream  energy  flow)  correspond¬ 
ing  to  the  center  of  each  cell.  The  net  momentum  and  energy 
transfer  to  the  body  is  also  stored,  and  from  this  the  drag 
coefficient  and  other  aerodynamic  quantities  can  be  calculated. 
The  contents  of  various  other  counters  showing  the  progress  of 
the  calculation  are  also  printed,  e.g.,  number  of  collisions, 
average  free  path  between  collisions,  etc. 

After  the  sampling  has  been  completed,  one  could  simply 
resume  the  calculation.  However,  if  this  is  done,  observations 
on  the  test  particle  taken  long  ago  will  be  weighted  equally 
with  the  latest  observations,  even  though  the  latest  are  pre¬ 
sumably  based  on  an  improved  estimate  of  the  distribution  func¬ 
tion  of  the  scatterers.  On  the  other  hand,  the  cumulative 
registers  may  not  be  reset  to  zero,  since  then  all  information 
based  on  earlier  calculation  would  be  lost.  The  procedure 
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adopted  here  is  to  decide  (perhaps  somewhat  arbitrarily)  the 
relative  weights  to  be  attached  to  prior  and  new  data.  The 
cumulative  registers  in  every  cell  (but  not  the  numbers  of  the 
list  representing  individual  particles)  are  multiplied  by  some 
number  less  than  unity,  such  that  they  will  appear  to  have  been 
generated  by  a  smaller  number  of  particles  than  actually  was 
used.  The  details  of  this  "downgrading  of  prior  information" 
will  be  explained  later. 

Generally,  it  is  desired  to  perform  the  calculation  for  a 
series  of  Knudsen  numbers  at  a  fixed  stream  speed.  The  calcu¬ 
lations  are  performed  in  decreasing  order  of  Knudsen  number, 
always  starting  with  the  free  molecular  case  (infinite  Knudsen 
number).  At  the  start,  all  cells  are  empty,  and  there  is  thus 
no  estimate  of  the  distribution  function  available,  and  no 
collisions  can  be  calculated.  This  poses  no  difficulty  for  free 
molecular  flow,  since  the  particles  do  not  interact  and  the 
collision  calculations  can  be  bypassed.  When  the  free  molecular 
solution  has  been  found,  it  serves  as  a  first  approximation  to 
the  solution  for  the  highest  finite  Knudsen  number,  etc.  As  the 
Knudsen  number  is  lowered,  the  convergence  will  become  slower, 
since  more  collisions  will  be  calculated,  and  eventually  limita¬ 
tions  on  computer  time  will  force  a  halt  to  the  computation. 

The  program  is  written  to  compute  collisions  by  either  the 
hard-sphere  or  the  "Maxwellian"  (inverse  fifth  power  repulsion) 
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molecular  model.  A  subroutine  could  also  be  written  to  compute 
collisions  by  some  empirical  model,  if  the  two  simple  models 
are  considered  too  unrealistic. 

Provision  is  made  for  writing  all  variable  values  on  a 
tape,  so  that  the  computation  could  be  continued  at  a  later 
time  after  examination  of  the  data. 

The  mathematical  basis  of  the  calculation  is  presented  in 
Section  IV.  When  results  are  available,  they  will  be  presented 
in  a  subsequent  report. 
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IV.  MATHEMATICAL  BASIS  OF  THE  SIMULATION 

In  this  section  the  mathematical  basis  of  the  various 

parts  of  the  simulation  will  be  presented.  Much  of  the  same 

material  is  presented  somewhat  more  generally  elsewhere,  es- 

8  5  4 

pecially  by  Patterson  ,  Haviland  ,  and  Lavin  .  The  treatment 
here  is  intended  as  a  derivation  and  explanation  of  the  method 
used  in  the  program  and  thus  follows  it  closely  in  order  and 
in  nomenclature. 

(1)  Dynamics  of  elastic  collisions 

Consider  two  particles  of  equal  mass  and  velocities  v^ 
and  Vg  which  collide  elastically.  Their  velocities  after 
collision  are  v^ '  and  Vg * .  Their  relative  velocity  before 
collision  is 


vr  V1  v2  * 


(1) 


From  conservation  of  momentum, 


V1  +  v2  =  vl'  +  v2*’ 


(2) 


and  from  conservation  of  energy, 


2  2  ,2  ,2 

V1  +  v2  =  V1  +  v2  » 


(3) 


where  v^  =  | v^ | ,  etc.  Combining  equations  (2)  and  (3),  we  obtain 
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v  1  =  V 

r 


r 


v1»  =  l/2(v1  +  v2)  +  1/2  vr' 


(4) 


v2’  =  l/2(v1  +  v2)  -  1/2  vr*. 


Note  that  equations  (4)  are  seven  scalar  equations  in  nine  un¬ 
knowns.  The  remaining  two  unknowns  correspond  to  the  direction 
of  the  new  relative  velocity.  This  depends  on  the  particular 
collision  model,  on  the  collision  parameters  and  on  v  . 

(2)  Collision  mechanics  for  a_  general  repulsive 


force  law 


Assume  a  repulsive  force  K/rn  per  unit  mass  between  the 


two  molecules,  where  r  =  -  x2 .  Then  the  equation  of  motion 

is  easily  shown  to  be 


r  =  2Kr/rn+1 


(5) 


In  polar  coordinates  these  equations  may  be  written  (see 
Figure  2) 


,  /0,  .2  2.2x  2K  1-n  .  ,  /0  2 

l/2(r  +  r  9  )  +  r  =  const  =  1/2 


(6) 


r20  =  const  =  bvr, 
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Figure  2. 


The  coordinate  system  in  the  collision  plane 
for  a  repulsive  force  law. 
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where  b  is  the  distance  by  which  they  would  miss  each  other  if 
no  interaction  occurred. 

Eliminating  time  from  (6), 


dr 

d§ 


2 

r 


4Kr 


5-n 


(n-l)b2v  2 
r 


(7) 


Introducing  the  dimensionless  coordinate 


W  =  b/r, 

v  2  1/n-l 

and  letting  WQ  =  bC-^-) 
equation  (7)  takes  the  form 


dW 

de 


2 

n-1 


(JL) 

o 


n-1 


(8) 

(9) 


(10) 


Now  at  the  point  of  closest  approach 

dW  dW  dr  n 

d9  dr  de 

This  is  the  intersection  with  the  apse  line,  and  corresponds  to 
the  angle  <j>  in  Figure  2  and  to  W  =  W  .  Thus  we  can  integrate 
(10)  to  get 
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w  n-1  -1/2 

(f) 

o 


dw  =  <>(wo) , 


(ii) 


w 


2 


2 

n-1 


where  W^(Wq)  is  a  positive  root  of 


1 


W 


2 

n-1 


(JL) 

vw  ' 

o 


n-1 


0. 


(12) 


The  collision  angle  (the  angle  between  the  new  and  old  relative 
velocities)  is  then 


X(WQ)  =  tt  -  2<j)(Wo).  (13) 

(3)  Maxwellian  molecule 

When  n  =  5,  the  above  equations  assume  a  simple  form.  A 
molecule  with  this  force  law  is  called  a  "Maxwellian"  molecule. 
Note  that  Wq  is  a  function  of  the  magnitude  of  the  relative 
velocity  and  of  the  collision  parameter  b  (the  "miss  distance") . 

Setting  n  =  5,  equation  (11)  can  be  integrated.  The  colli¬ 
sion  angle  then  is 

x(WQ)  -  TT - -  -4-r 74  Kf-l/2  -  1/2(1  +  ~^)~1/21,  (14) 

°  (1+2/W  4)i/4  L  W  4  J 

o  o 

where  K(a)  is  the  complete  elliptic  integral  of  the  first  kind: 
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1  3-35-103191 


Figure  3.  Collision  geometry  for  hard  sphere  molecules. 
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(15) 


tt/2 


K(a)  -  J 


.  2 


-  asm 


y 


(4)  Hard  sphere  molecule 

The  "hard  sphere"  molecule  is  a  limiting  case  of  the  in¬ 
verse  power  molecule.  The  scattering  angle  x  then  becomes  inde¬ 
pendent  of  the  speeds  of  the  particles,  depending  only  on  the 
"miss  distance"  b.  It  is  easily  seen  from  the  diagram  of  the 
collision  in  Figure  3  that 


b  =  R  cos (x/2) 


(16) 


(5)  Calculation  of  velocities  after  collision 
For  both  Maxwellian  and  hard  sphere  molecules,  it  has  been 
shown  above  how  to  calculate  X,  the  angle  between  the  new  and 
old  relative  velocities,  given  the  magnitude  of  the  relative 
velocity  vr  and  the  "miss  distance”  b  (only  the  latter  is  re¬ 
quired  for  the  hard  sphere  molecule).  In  addition  to  the  scatter¬ 
ing  angle  X,  there  is  another  angle,  e,  which  must  be  known  before 
the  collision  is  completely  specified.  It  is  the  angle  between 
the  plane  containing  the  new  and  old  relative  velocities  and  some 
arbitrary  plane  containing  the  old  relative  velocity.  The  range 
of  x  Is  from  0  to  rr,  while  c  varies  from  0  to  2tt .  The  calcula¬ 
tion  of  e  is  very  simple  for  any  central  force  law.  With  any 
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isotropic  scattering  model,  all  values  of  e  from  0  to  2tt  are 
equally  likely. 

Once  e  and  X  are  known,  the  new  velocities  can  be  computed 
as  follows: 

Form  a  set  of  three  unit  vectors  z^,  z,^,  and  z^  such  that 
z^  is  parallel  to  the  old  relative  velocity.  The  new  relative 
velocity  is  then 


v  •  =  v  (z,  cosX  +  z0  sinX  cose  +  z„  sinX  sine)  (17) 

r  i*  J.  z  u 


The  unit  vectors  z^ ,  z^,  and  z^  may  themselves  be  expressed  in 

terms  of  z  ,  z  ,  and  z  ,  and  equation  (4)  may  then  be  used  to 
x  y  z 

compute  the  new  velocities  v^ '  and  ' .  This  applies  to  any 
scattering  model . 

An  even  easier  procedure  is  available  for  the  hard  sphere 
molecule,  based  on  the  peculiar  fact  that  in  a  hard  sphere 
collision  the  direction  of  v  '  is  independent  of  the  direction 
of  v  ,  and  hence  may  be  chosen  randomly. 


(6)  Selection  of  collision  parameters 

Suppose  b  „  is  the  maximum  miss  distance  for  which  colli- 
sions  are  to  be  calculated.  Then,  in  the  plane  containing  the 
target  particle  and  normal  to  the  trajectory  of  the  incident 
particle,  let  a  circle  of  radius  2b  centered  at  the  target 
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Figure  4.  Geometry  of  collision  parameters  for  binary 
collisions.  Orbit  lies  in  plane  containing 
target  particle. 
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particle  be  drawn  (Figure  4) .  If  it  has  been  decided  that  the 
two  particles  shall  collide,  then  the  trajectory  of  the  incident 
particle  (in  the  coordinate  system  of  the  target  particle)  must 
intersect  this  infinitesimal  circle.  Clearly  there  can  be  no 
variation  in  the  distribution  function  over  such  a  small  circle, 
so  all  intersection  points  are  equally  likely.  If  and  R2 
are  two  rectangular*  random  numbers,  a  random  point  inside  the 
circle  is  given  by 


e  =  2ttR-^  , 


b2  =  b2  R„. 
max  2 


(18) 


This  is  the  procedure  for  choosing  collision  parameters. 


( 7 )  Calculation  of  collision  rates 

In  the  Monte  Carlo  scheme,  a  test  molecule  enters  a  cell, 
and  we  must  decide 

(1)  whether  or  not  a  collision  will  occur,  and 

(2)  in  the  event  of  a  collision,  the  velocity  of 
the  collision  partner. 

To  make  these  decisions  we  must  know  the  scattering  model,  the 

"Rectangular"  random  numbers  are  uniformly  distributed  between 
0  and  1. 
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velocity  of  the  test  molecule,  and  the  density  and  velocity 
distribution  of  "target"  molecules  in  the  cell. 

Now  consider  two  beams  of  particles  of  velocities  v^  and 
Vg  which  intersect  such  that  their  collision  parameters  are 
within  the  element  area  b  db  de,  in  the  coordinate  system  of 
the  target  beam.  Then  in  a  time  interval  dt ,  all  the  particles 
in  the  incident  beam  lying  within  the  volume  element  vr  b  db  dt  de 
will  collide.  Summing  over  both  collision  parameters  gives  us 
the  collision  integral: 


X(  v 


1’ 


d  e 


J^max 
J  o 


v  b  db 
r 


(19) 


This  is  the  volume  per  unit  time  of  the  incident  beam  which 
collides  with  the  target  beam. 

For  the  hard  sphere,  b  is  obviously  twice  the  radius  of 
the  sphere.  For  other  collision  models,  the  choice  is  not  ob¬ 
vious,  since  collisions  can  occur  for  any  value  of  b.  Normally 
we  do  not  wish  to  waste  time  computing  grazing  collisions,  so 
we  choose  b  „  such  that  collisions  with  deflection  angle  X  less 
than  some  minimum  value  (usually  about  11°)  will  be  ignored. 

This  generally  makes  b  a  function  of  velocity.  For  the 

max 

Maxwellian  molecule,  equation  (14)  shows  that  it  is  W  which 

o 

has  a  fixed  maximum  value  (independent  of  velocity)  for  a  fixed 

minimum  scatter  angle.  Thus  instead  of  choosing  b  =  b  R, ,  as 

max  1’ 
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in  (18),  it  is  more  convenient  to  choose  W  =  W  „  R,  . 

o  o  max  i 

The  collision  rate  a(v^)  is  then  obtained  by  integrating 
over  the  distribution  f(v2)  of  target  particles: 

a(vi>  =  J  dv2  fCvgJXCv^^,  v2).  (20) 

a11  ^2 

For  hard  sphere  molecules,  the  collision  rate  is  easily 
found  from  equations  (19)  and  (20): 

^(v^)  =  nr2  J  vrf(v2,x)  dv2 ,  (21) 

all  vx 

where  r  is  the  radius  of  the  hard  sphere. 

For  the  Maxwellian  molecule,  we  use  equations  (9)  (with 
n  =  5)  and  equations  (19)  and  20)  to  obtain 

°<72>  -  P"»o2.a,  <2K)1/2’  <22) 

where  p(x)  =  f(v2,x)dv2  is  the  local  number  density, 

‘’all  v^ 

5 

It  was  shown  by  Haviland  that  a  value  of  Wq  max  of  about  1.5 
gave  good  results  (this  corresponded  to  a  cut-off  of  about  11° 
for  the  scatter  angle) . 
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(8)  Mean  free  path 

In  the  Monte  Carlo  simulation,  the  stream  speed  and  the 
Knudsen  number  are  specified,  and  parameters  such  as  K  in 
equation  (22)  and  r  in  equation  (21)  are  not  meaningful.  To 
eliminate  these  parameters  from  our  equations,  we  consider  the 
mean  free  path  X. 

g 

Patterson  shows  that  the  mean  free  path  for  hard  spheres 
of  radius  r  is 


X 


_ 1 _ 

%/2nr2p 


(23) 


while  for  Maxwellian  molecules, 


X 


2 


3/2 

TT 


pW 


2 

o  max 


1/2 


(24) 


Now,  since  WQ  max  is  arbitrary,  this  gives  an  arbitrary 

5 

mean  free  path.  Haviland  shows  that  the  "effective  mean  free 
path"  (based  on  comparisons  of  viscosity  and  heat  conduction 
with  other  collision  models)  is 


'ef  f  2tt(0.436)p 


(— ) 
v  K  ' 


1/2 


(25) 


Thus,  if  we  want  to  specify  a  certain  mean  free  path  (for  com¬ 
parison  with  other  models)  X^^,  we  choose 
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kT 


1/2 


(-£-)  =  2n(0 . 436)  pXe , 


(26) 


The  mean  free  path  in  the  computation  process  then  is  (for  the 
Maxwellian  molecule) 


4(0.436)X 


X  = 


eff 


(27) 


ttW 


o  max 


(9)  Representation  of  the  distribution  function  of  a^ 
finite  collection  of  particles 

In  the  present  Monte  Carlo  scheme  the  velocity  distribution 
is  represented  by  a  set  of  N  "particles"  with  velocities  and 
associated  weighting  factors  w^ .  At  a  given  instant,  this  is  a 
crude  representation,  but  in  the  long  run  an  ensemble  average 
of  these  collections  should  give  a  good  approximation  to  the 
true  distribution  function. 

Assuming  a  collection  of  particles  has  been  "correctly" 
selected,  we  focus  our  attention  on  computation  of  the  collision 
rate,  and  when  necessary,  the  selection  of  a  collision  partner. 

At  any  instant,  our  estimate  of  the  distribution  function 
is  a  finite  collection  of  "spikes": 


Wv)  ■ 


N 

I 


est  ;  i 
i^i 


w.  6 ( v  -  v.) 


(28) 


y  w. 


i=l 
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where  6(v)  is  the  usual  Dirac  delta  function,  i.e., 


JIT  dvxdvydV(7)  ,  !. 

all  v 


Using  (28)  the  collision  rate  for  hard  sphere  molecules, 
when  the  test  particle  has  velocity  v,  is 


a(  v) 


N 

2  i=l 
TTr  pest  N 


i=l 


(29) 


Far  away  from  the  body  (free  stream  conditions)  we  can  write 
(23)  as 


r 


2 


f 


where  X  is  the  free 
written 


a  ( v) 


pest 

PQ 


stream  mean  free  path,  and  (2)  can  then  be 


1 

\s/2 


N 

y>.iv  - 

i=l 


N 

y  w. 

L  x 
i=l 


(30) 


This  is  the  collision  rate  for  hard-sphere  molecules.  If  Kn  is 
the  free  stream  Knudsen  number  and  d  is  the  body  diameter,  then 
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X  =  Knd 


and  all  that  is  lacking  is  a  means  for  estimating  the  ratio  of 
the  local  density  to  the  free-stream  density.  This  last  point 
will  be  taken  up  later.  (See  equation  (40).) 

The  collision  rate  for  Maxwellian  molecules  is  independent 
of  all  velocities  (this  is  the  motivation  for  the  Maxwellian 
collision  model)  and  is  obtained  by  combining  equations  (22) , 
(24)  and  (27) : 


0=0  ,  nW2  (2K)1/2  M)1/2 

est  o  max  kT 


(31) 


w 


o  max 


2(0. 436) X 


eff 


V2kT 


One  of  the  normalizations  used  in  the  Monte  Carlo  calcula¬ 
tion  is  to  set  the  mean  thermal  speed  (the  root-mean  square 
deviation  from  the  mean  speed)  equal  to  3/2.  This  is  equivalent 
to  setting 


kT  =  1 


(32) 


and  when  (32)  is  substituted  into  (30),  the  collision  rate  (per 
"incident"  molecule)  for  Maxwellian  molecules  becomes 
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W  o 

o  max 


(33) 


(0.436)Xeff  y/2 


It  should  be  remembered  that  the  "mean  free  path"  is  not  a  well- 
defined  quantity  for  all  collision  models.  It  is  readily  vis¬ 
ualized  for  the  hard-sphere  models,  but  for  any  inverse  power 
law  molecule  (such  as  the  inverse  5th  power  "Maxwellian")  the 
definition  is  quite  arbitrary.  For  the  Maxwellian  molecule, 
the  effective  mean  free  path,  defined  in  (25)  is  used  exclusively 
here.  With  this  choice,  for  a  given  Knudsen  number  and  with  all 
other  parameters  identical,  the  hard  sphere  and  Maxwellian  models 
will  have  the  same  heat  conduction  and  viscosity  in  the  limit  of 

continuum  flow.  One  must  be  very  careful  not  to  attach  any  other 
significance  to  the  mean  free  path  or  to  the  "Knudsen  number" 
when  comparing  calculations  based  on  different  collision  models. 

(10)  Deciding  when  the  next  collision  will  occur 

For  any  binary  collision  model,  the  probability  that  a 
given  molecule  of  velocity  v  does  not  suffer  a  collision  in  time 
interval  t  is  exp  [-a(v)t]  where  ct(v)  is  the  local  collision 
rate,  assumed  constant  over  the  small  distance  vt .  The  proba¬ 
bility  density  function  for  the  next  collision  occurring  at  t 

c 


is  then 


p(t  )  =  ct  exp (-at  ). 


(34) 
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We  now  (with  some  foresight)  make  the  transformation 


t  =  -  -  In  R. 

C  CT 


(35) 


The  probability  density  function  for  R  is  then  easily  seen  to 
be 


q(R)  =  p(tc)|-5|-|  =  1.  (36) 

Since  the  range  of  t  is  0  <  t  <  “,  the  range  of  R  must  be 

c  c 

0  <  R  <  1,  and  thus  R  is  a  "rectangular  random  number".  Thus 
we  can  draw  a  collision  time  tc  by  simply  drawing  a  rectangular 
random  number  and  using  (35). 

(11)  Choice  of  collision  partner 

If  the  collision  time  t  turns  out  to  be  less  than  the 

c 

time  t  required  to  traverse  a  given  cell,  then  a  collision  will 
w 

occur  in  that  cell.  It  is  then  necessary  to  select  a  collision 
partner,  so  that  the  new  velocity  of  the  test  molecule  may  be 
computed.  Given  the  collision  model,  the  local  velocity  distri¬ 
bution  and  the  velocity  of  the  incident  particle,  one  can  always 
write  the  distribution  function  for  the  velocity  of  the  colli¬ 
sion  partner,  conditional  on  the  fact  that  a  collision  does 
occur.  One  then  draws  a  collision  partner  from  this  distribution. 
(This  does  not  in  any  way  alter  the  distribution  function.) 
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Now  the  probability  of  the  test  particle  colliding  with 


a  particle  of  velocity  (given  that  a  collision  occurs)  is 
clearly  proportional  to  the  contribution  of  the  velocity  v^  to 
the  total  collision  rate.  Thus  we  arrive  at  a  simple  rule  for 
selecting  a  collision  partner: 


(a)  choose  a  rectangular  random  number  Rn; 

(b)  choose  the  velocity  v’k  which  satisfies  the  inequality 

k-1  N  k 

V  W.|v  -  v.|<Rn  Y  wi|v  -  vjs  Y  w.|v  -  v.|  (37) 

i=l  i=l  i=l 

for  the  hard-sphere  molecule;  or 

(c)  for  the  Maxwellian  molecule, 

k-1  N  k 

w.  ^  Y  w.  .  (38) 

1  L  i 

i=l  i=l  i=l 


Y  w.  <  r  Y 

L  i  n 


By  this  procedure,  the  test  particle  will  be  scattered  in  the 
same  way  as  would  a  particle  in  the  real  gas,  provided  that  the 
delta  function  representation  of  the  velocity  distribution  is 
an  unbiased  sample  from  the  true  distribution. 


(12)  Estimation  of  the  distribution  function 

So  far  we  have  calculated  the  behavior  of  the  test  particle, 
assuming  that  the  distribution  function  is  correctly  represented 
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by  a  set  of  weighted  delta  functions  in  velocity-space,  and 
that  some  reasonable  estimate  of  the  local  number  density  is 
available . 

We  now  reverse  this  view  and  assume  that  the  test  particle's 
behavior  is  representative  of  particles  in  the  actual  gas,  i.e., 
that  the  true  distribution  function  can  be  deduced  by  observing 
the  test  particle  for  an  infinitely  long  time. 

Now,  each  cell  in  the  x-y  grid  can  be  regarded  as  a  window 
which  occasionally  observes  the  test  particle.  If  on  the  ith 
passage  through  the  cell  in  question  the  test  particle  has  veloc¬ 
ity  v.  and  spends  time  w^  in  traversing  the  cell,  then  clearly 
the  best  estimate  (given  no  other  information)  that  the  cell  can 
make  of  the  velocity  distribution  is  a  sum  of  delta  functions  at 
positions  v\  with  areas  w^.  This  is  simply  all  the  cell  "sees"; 
any  other  representation  would  involve  some  prior  notion  of  the 
form  of  the  distribution  function. 

The  estimation  of  the  ratio  of  local  number  density  to 
free-stream  number  density  (hereafter  called  normalized  density) 
must  also  be  based  on  observations  of  the  test  particle.  For 
this  estimate  we  take 

est  _  (time  test  particle  has  spent  in  given  cell)  (39) 

0  (time  it  would  have  spent  in  same  cell  under 

0  free  stream  conditions) 
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The  numerator  in  (39)  is  just  the  sum  of  all  w^ .  The  denom¬ 
inator  is  not  so  obvious.  One  can  not  use  the  total  test  par¬ 
ticle  time,  as  this  would  force  the  average  density  over  the 
x-y  grid  to  be  the  free-stream  density,  which  is  obviously  in¬ 
correct.  Instead,  we  count  the  number  of  particles  N  ^  which 
have  entered  at  the  ends  (upstream  or  downstream)  and  multiply 
by  the  expected  time  that  a  random  upstream  particle  should 
spend  in  the  given  cell.  Letting  Ix,  I  be  the  integer  coor¬ 
dinates  designating  the  given  cell  in  the  grid,  our  estimate  of 
the  number  density  is  then 


N 

(  Y  w .  )  L  L 

L  1  x  y 

pest  _ 

p  N  .  ACE  ,1  ) 

Ho  pi  x’  y 


U_ 

L 


(40) 


where  U  is  the  stream  speed,  L  and  L  are  the  over-all  dimen- 

x  y 

sions  of  the  grid,  and  A(Ix,I^)  is  the  area  of  cell  (Ix,I^t). 

It  was  mentioned  in  Part  III  that  the  list  of  delta  func¬ 
tions  in  a  given  cell  was  never  allowed  to  grow  beyond  some 
maximum  length,  (hereafter  called  Ig)  .  However,  Y  w^  is  also 
accumulated  in  a  register  for  each  cell,  so  that  density  is 

estimated  from  all  the  times  the  test  particle  has  entered  the 
cell  (since  the  last  time  the  register  was  reset  to  zero),  not 
just  from  the  current  members  of  the  list. 
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Similarly,  the  moments  of  the  distribution  function  which 
are  of  special  interest,  namely  the  velocity  and  energy,  are 
accumulated  in  registers.  These  are  not  involved  in  the  calcu¬ 
lation  of  collisions  but  are  printed  out  whenever  a  "sample''  is 
taken. 

The  "normalized  velocity  array",  which  is  printed  in  tabu¬ 
lar  form,  is  simply 

N 

Y  w.  v . 

L  11 

i=l 

V(I  ,1  )  =  -rr -  (41) 

x’  y  N 

I  *i 

i=l 

where  N  represents  all  the  particles  which  have  entered  the  cell 
(Ix,Iy),not  just  the  number  on  the  current  list.  The  normalized 
temperature  is 


N 


I  l7il  % 


T(  I  ,1  )  = 

x  y 


i=l 


N 


y 

i—i 

i=l 


-  |v(I  ,1  )| 
i  x’  y  1 


w. 

L,  1 


(42) 


The  normalized  stagnation  temperature  is 


N 


Tc ( I  ,1  )  = 
S  x’  y 


|l  I  Vi 

i=l 


N 


!>  i  ,t 

i=l 


(43) 


39 


It  is  also  possible  to  accumulate  other  moments,  if  they  should 
be  of  interest. 

(13)  Deletion  of  entries  from  the  lists  of  delta  functions 

Every  time  the  test  particle  traverses  the  cell,  informa¬ 
tion  is  gained  about  the  local  distribution  function.  This  is 
added  to  the  prior  information  as  described  in  Section  (12). 

When  the  list  of  delta  functions  has  reached  maximum  length 

I  ,  an  entry  must  be  deleted  to  make  room  for  each  new  entry. 
Each  entry  spends  the  same  time  on  the  list,  regardless  of  its 
weighting  factor.* 

(14)  Boundary  conditions 

When  the  test  particle  crosses  the  symmetry  plane  x  =  0, 
it  is  reflected  specularly,  i.e.,  the  signs  of  y  and  vy  are 
changed. 

When  it  crosses  the  end  or  side  boundaries,  it  is  lost  and 
a  new  test  particle  is  generated  by  the  method  described  in 
Section  ( 15) . 

When  it  strikes  the  body,  it  is  absorbed,  and  a  new  particle 
is  emitted  at  the  same  point  by  the  procedure  described  in  Sec¬ 
tion  ( 16)  . 


See  the  alternative  procedure  described  in  the  Appendix. 
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(15)  Generation  of  new  test  particles 
a.  generation  of  position 

When  a  new  test  particle  is  to  be  introduced,  its  position 
is  chosen  at  random  along  the  wide  or  end  boundary.  A  certain 
fraction  3  must  be  introduced  along  the  side  boundary  to  com¬ 
pensate  for  outward  diffusion.  Otherwise,  even  without  the  body, 
the  density  would  decrease  in  going  from  the  upstream  end  to  the 
downstream  end.  It  is  clear  that  the  particles  introduced  at 
the  sides  must  compensate  only  for  the  free-stream  diffusion, 
not  for  particles  deflected  sidewards  by  the  body.  The  side 
boundary  is  not  a  rigid  wall.  The  only  approximation  caused  by 
the  nearness  of  the  side  boundary  is  the  neglect  of  collisions 
outside  this  boundary.  In  free-molecular  flow,  the  boundary  has 
no  effect  at  all. 

The  fraction  3  is  calculated  as  follows: 


J 


s 


flow  thru  side 


2 


2  2i 

+  Vy  +  V2  ] 


(44) 


J  =  flow  thru 
e 


all  v 
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Integrating  equation  (44), 


J  = 
s 


1 


L 

x 


(45) 


J  =  L  {  —  exp[-u2]  +  u  erf (u) } 

y  A 

Thus,  the  fraction  of  particles  entering  at  the  sides  is 


P  = 


J  +  J 
s  e 


L  L 

1  +  2  exp(-u^)  +  2 /rr  uerf(u) 

X  x 


(46) 


3  ^ 


1  +  V"  L*  U 

x 


if  u  s  3. 


The  position  of  the  new  test  particle  is  then  determined  as 
follows : * 

(1)  draw  R  ; 

n  R 

(2)  if  R  <0,  introduce  at  side,  a  distance  — £■  L  from 

n  ’  3  x 

the  upstream  end; 

(3)  if  Rn  >  3j  introduce  at  upstream  end  if  velocity  is 
positive,  downstream  if  velocity  is  negative,  at 

j|( 

This  procedure  will  be  replaced  by  a  stratified  sampling  pro¬ 
cedure,  designed  to  reduce  the  effect  of  uneven  flow  of  particles 
across  the  upstream  end. 
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3 


y  = 


R  - 
n 

1-3 


L  . 

y 


The  new  velocity  must  be  drawn  from  the  free  stream  veloc¬ 
ity  distribution: 


f(v)  =  (^)3//2  exp{-[  ( v  -  u) 

A 


2  2  2,, 
+  v  +  v  1 } . 
y  z  J  J 


(47) 


To  choose  a  normal  random  deviate,  we  use  a  very  efficient  method 

9 

suggested  by  Kahn: 

(1)  choose  random  numbers  R^ ,  R^', 

(2)  let  A  =  -lnR^,B  =  -In  Rg ; 

(3)  if  (A  -  l)2  >  2B,  go  back  to  1; 

2 

(4)  if  (A  -  1)  s  2B,  accept  A,  which  is  then  a  normal 
deviate,  with  variance  1. 

(5)  Choose  at  random  the  sign  to  be  affixed  to  A. 

To  show  that  this  procedure  does  indeed  yield  a  normal  random 
variate,  we  compute  the  probability  density  function  for  A, 
given  that  A  is  accepted.  By  Bayes*  rule, 


P[A | (A  -  l)2  S  2B]  -  P[(A  -  l)2  «  2B|AJP(A) 


(48) 


P[ (A  -  1)"  s  2B] 


where  P(C|D)  denotes  the  "probability  of  event  C,  given  event 
D",  and  P(C|D)  denotes  the  "probability  density  function  for 
quantity  C,  given  event  D". 
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The  denominator  of  the  right  side  of  (48)  need  not  be 
evaluated  as  it  is  not  a  function  of  A,  but  only  a  constant 
which  normalizes  the  density  function.  From  Section  (10)  and 
from  step  (2)  above,  we  recognize  that  A  and  B  have  exponential 
distributions : 


P(B)  =  exp  (-B),  P(A)  =  exp  (-A). 


(49) 


Finally  the  remaining  factor  in  (48)  is  easily  evaluated  as  a 
function  of  A: 


2 


exp  ( -B) dB 


(A-l) 


2 


(50) 


Combining  (48),  (49)  and  (50): 


p[A|(A  -  l)2  £  2B]  =  C1  exp  f- 


(A-l) 

2 


2 


■]  exp  ( -A) 


=  C2  exp  (-  ^-) 


(51) 
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This  is  the  desired  normal  distribution  for  A,  and  C ^  is  the 
usual  normalization  constant.  This  method  for  generating  normal 
variates  is  76  percent  efficient  (i.e.,  it  rejects  only  24  per¬ 
cent  of  the  pairs  of  rectangular  random  numbers)  and  is  probably 
the  fastest  method  for  most  digital  computers. 

Since  we  want  the  variance  of  each  of  the  three  components 
to  be  1/2,  we  multiply  the  numbers  obtained  by  Kahn's  method  by 
1A/2.  The  stream  speed  U  is  then  added  to  the  x-component  of 
velocity . 

The  number  of  particles  introduced  at  the  ends  and  the 
sides  (N  ^  and  Nside>  respectively)  are  recorded  by  a  counter. 

N  .  is  used  to  determine  when  to  print  out  and  when  to  end  the 

pi 

computation.  Other  counters  record  the  energy  and  momentum  of 
the  newly  generated  particles. 


(16)  Interaction  of  particles  with  the  aerodynamic  surface 
This  section  describes  the  way  in  which  the  test  particle 
is  re-emitted  when  it  strikes  the  body.  In  the  present  calcula¬ 
tion  the  body  is  a  right  circular  cylinder,  but  any  body  with 
cylindrical  symmetry  (i.e.,  no  variation  in  the  z-direction) 
could  be  treated  by  this  method.  Similarly,  other  surface  inter¬ 
action  models  could  easily  be  substituted  for  the  present  model. 

In  the  present  calculation,  a  fraction  of  the  particles 
striking  the  surface  are  reflected  specularly  (i.e.,  the  compo¬ 
nent  of  velocity  normal  to  the  surface  is  reversed) ,  and  the 
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rest  are  replaced  by  new  particles  emitted  diffusely  at  a  tem¬ 


perature  such  that  the  thermal  speed  is  times  the  free  stream 
thermal  speed. 

One  may  conceptually  visualize  a  small  hole  in  the  surface 
at  the  impact  point,  behind  which  is  a  gas  in  equilibrium  with 
thermal  speed  S^.  A  stream  of  particles  diffuses  outward  from 
this  hole.  We  then  select  at  random  one  of  these  particles  to 
be  the  new  test  particle.  Clearly  the  probability  density 
governing  this  selection  must  be  proportional  to  the  flux  of 
particles  of  a  given  velocity  through  the  hole.  This  flux  is 
proportional  to  the  velocity  distribution  times  the  velocity 
component  normal  to  the  hole. 

Let  the  normal  and  the  two  tangential  components  of  veloc¬ 
ity  be  vn,  v.^,  vt2>  respectively.  The  density  function  for 
particles  passing  through  the  hole  is  then 


P(v 


n 


Vtl’Vt2)  ' 


lv 
1  n 


exp 


+  v 


yi 


+  V 


t22]} 


(52) 


where  An  is  the  normalization  constant,  and  v  >0. 

1  n 

The  normal  component  is  then  selected  according  to  the 
following  procedure: 


P(v  ) 
n 


A2  exp 
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2 


v 


Let  u 


n 


S 


2  • 


b 


Then  the  probability  density  function  for  u  is 


q(u)  -  p[vn(u)-j£] 


=  Ag  exp( -u) . 


But  we  have  already  seen  in  Section  (10)  how  to  draw  from  this 
distribution,  given  a  rectangular  random  number  R^.  We  simply 
take 

u  =  -lnR^ 

and  the  normal  component  of  velocity  then  becomes 


(53) 


The  two  tangential  velocity  components  v^^  and  v ^  could  be 
chosen  by  the  method  of  Section  (15).  However,  it  is  more  con¬ 
venient  to  choose  the  magnitude  of  the  tangential  velocity  vector, 
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and  then,  recognizing  that  all  angles  in  the  tangential  plane 
are  equally  likely,  choose  a  random  angle  between  0  and  2rr. 
Expressed  in  these  polar  coordinates,  the  probability  density 
for  is  then 

2 

P(vt)  =  A4vt  exp  [ - ^2~\  (54) 

Sb 

i.e.,  v.  and  v  have  the  same  density  functions.  Thus  we  derive 
t  n 

— ♦ 

vt  from  two  rectangular  random  numbers: 


vt  *  Sb'/'  lnR2 

<{)  =  2TrRg  . 


(55) 


where  it  is  convenient  to  measure  from  a  line  parallel  to  the 
cylinder  axis.  The  new  velocity  components  in  x-y-z  coordinates 
are  then  obtained  by  a  straightforward  transformation. 


( 17)  Sampling 

After  a  pre-assigned  number  Z^sam  of  test  particles  has 
been  generated  at  the  end  boundaries,  a  sample  is  taken. 

The  current  values  of  the  arrays  described  in  Section  (12) 
equations  (39)  -  (43)  are  printed  in  such  a  way  that  the  flow 
field  may  be  visualized.  The  contents  of  the  various  counters 
are  also  printed. 
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Various  aerodynamic  coefficients  are  also  calculated  and 
printed  out: 

The  drag  coefficient  CD  is  the  total  momentum  per  cross- 
section  area  transferred  to  the  cylinder,  divided  by  the  total 
momentum  per  unit  area  which  has  entered  upstream  since  the 
previous  sampling.  A  factor  2  is  added  to  conform  to  the  con¬ 
vention  that  Cp  =  2  if  all  the  stream  momentum  over  the  cross 
section  of  the  body  is  absorbed  by  the  body. 

The  energy  spread  of  the  entering  particles  is  calculated, 
mainly  as  a  check  on  the  generating  process. 

The  average  free  path  for  all  test  molecules  since  the 
previous  sampling  is  also  calculated.  This  may  be  expected  to 
be  somewhat  lower  than  the  free-stream  mean  free  path,  due  to 
increased  density  in  front  of  the  body. 

The  energy  transfer  to  the  cylinder,  normalized  by  the 
energy  that  would  be  absorbed  under  free  stream  conditions  if 
all  particles  within  the  cross  section  of  the  cylinder  were 
absorbed . 

Finally  the  three  components  of  force,  per  entering  par¬ 
ticle,  exerted  on  the  cylinder  are  calculated,  mainly  as  a  check 
on  statistical  fluctuations. 

After  the  sampling  is  completed,  the  calculation  may  be 
continued  at  the  same  or  at  a  lower  Knudsen  number,  until  the 
number  of  particles  entering  upstream  Z  ^  exceeds  the  next 
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pre-assigned  number  Z  .  .  However,  if  this  process  were  con- 

tinued,  the  quantities  which  have  been  accumulated  in  each  cell 
would  soon  become  insensitive  to  changes  in  Knudsen  number, 
since  the  amount  of  information  derived  from  a  single  passage 
of  the  test  particle  would  become  very  small  in  comparison  with 
the  information  already  stored  in  the  cell.  To  remedy  this 
situation,  the  following  procedure  is  followed  immediately  after 
completing  a  sampling: 

(1)  One  "pretends"  that  the  latest  sample  (which 

was  actually  based  on  Z  .  entering  particles) 

J  pisam 

was  based  on  only  zpieff  particles,  where  zpieff  < 

Z  .  .  (Both  these  numbers  are  from  lists  sup- 

pxsam 

plied  on  data  cards.) 


(2)  Then  Z  .  and  all  the  current  accumulated  values 

pi  - 

in  the  array  of  cells  are  multiplied  by  zpieff/ 

Z  . 
pisam 

(3)  The  new  (lower)  Z  ^  is  now  considered  to  be  the 
number  of  particles  upon  which  all  present  esti¬ 
mates  are  based. 


(4)  The  calculation  then  proceeds  in  the  normal  way. 
Note  that  by  this  procedure  we  can  attach  any  desired 
weight  to  information  gained  prior  to  the  latest  sample. 
The  lists  of  individual  delta-functions  ("particles") 
in  each  cell  are  not  affected  by  this  procedure. 
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V.  CONCLUSIONS 


Most  of  the  main  features  of  this  simulation  have  been 
described  in  Part  IV.  Results  are  not  yet  available;  they  will 
be  presented  in  subsequent  reports.  It  is  also  probable  that 
other  features  will  be  added  to  the  program  to  make  the  model 
gas  more  realistic  and  to  consider  other  bodies  with  cylindrical 
symmetry,  such  as  a  flat  plate.  It  is  also  planned  to  perform 
a  similar  simulation  for  problems  with  axial  symmetry,  such  as 
a  sphere  or  a  cone.  All  such  refinements  will  be  described  in 
subsequent  reports. 
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APPENDIX 


An  Alternate  Representation  of  the  Velocity  Distribution 

In  the  model  described  above,  the  velocity  distribution  was 
represented  in  each  cell  by  a  fixed  number  of  weighted  delta 
functions : 

I 

max 

f ( v )  =  A  V  w  6(v  -  v  )  (A-l) 

j  JL  l 

i=l 

The  weighting  factors  w^  were  necessary  to  take  proper  account 
that  different  entries  v\  represent  different  intervals  of  time 
during  which  the  test  particle  was  observed.  By  making  w^  pro¬ 
portional  to  the  time  the  test  oarticle  spent  in  the  cell  during 
the  associated  passage,  the  contribution  of  each  v^  to  the  cur¬ 
rent  estimate  is  properly  accounted  for.  Each  entry  v^  remains 

on  the  list  for  I  passages  of  the  test  particle,  after  which 

max 

it  is  deleted  to  make  room  for  a  new  entry.  We  shall  refer  to 
this  scheme  as  Version  I  in  the  following  discussion. 

An  alternate  representation  of  the  velocity  distribution  is 
to  let  all  the  w^  in  A-l  be  unity ,  so  that  all  the  delta  func¬ 
tions  have  equal  area.  If  this  is  done,  the  length  of  time 
(number  of  passages  of  the  test  particle)  an  entry  stays  on  the 
list  must  be  proportional  to  the  time  the  test  particle  took  to 
traverse  the  cell.  We  refer  to  this  method  as  Version  II. 


53 


13-35-10321] 


ACTUAL  CONTINUOUS 
DISTRIBUTION 


f  Iv) 


yv) 


VERSION  I 


Uv) 


VERSION  n 


PASSAGE  1 


PASSAGE  9 


PASSAGE  17 


PASSAGE  25 


PASSAGE  33 


HISTOGRAM 
FROM  ABOVE 
5  SAMPLES 


4it  h,  n 

H 

\  i  i  i  L  k  k  Tl 

1 

t  i  »  i  .  t  li 

III 

\  *  t  *  *  i  *  i 

mi  i 

fm 

•••*% 

•  . 

Figure  5.  Schematic  comparison  of  two  representations  of 

the  velocity  distribution.  (The  histograms  were 
not  actually  calculated.)  The  small  circles 
indicate  spikes  which  were  also  present  at  the 
last  sampling. 
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It  is  clear  that  in  both  versions,  the  contribution  made  by 
entry  v^  to  the  "total  representation"  is  proportional  to  the 
weighting  factor  w^  times  the  length  of  time  the  entry  is  on 
the  list.  (By  "total  representation"  we  mean  the  approximation 
to  the  true  velocity  distribution  obtained  by  averaging  the 
lists  over  a  long  time  interval.)  Thus,  in  the  long  run,  both 
versions  are  statistically  correct  representations  of  f(v). 
However,  at  a  given  instant,  the  lists  in  the  two  versions  will 
have  rather  different  appearances.  To  illustrate  this  difference, 
Figure  5  shows  a  hypothetical  one-dimensional  velocity  distribu¬ 
tion  with  two  peaks,  and  possible  appearances  of  the  lists  (con¬ 
sisting  of  eight  entries)  at  successive  passages  of  the  test 
particle  are  shown  for  the  two  versions.  The  "histograms"  would 
be  calculated  by  dividing  the  v-axis  into  small  intervals  and 
taking  the  average  of  all  spikes  which  fall  within  an  interval. 
Note  that  in  both  versions  the  high  velocity  peak  is  represented 
in  better  detail  than  the  low  velocity  peak.  This  is  simply  be¬ 
cause  most  of  the  oncoming  particles  belong  to  the  high-velocity 
peak,  so  more  information  is  obtained  in  a  given  time  about  that 
part  of  the  distribution. 

It  would  appear  at  a  given  instant  that  Version  II  better 
represents  the  slow  particles  (i.e.,  treats  all  particles 
equally)  while  Version  I  favors  fast  particles,  giving  a  very 
crude  representation  of  the  low-velocity  peak.  While  this  is 
true,  it  is  misleading.  In  Version  II,  the  "slow"  entries  must 
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stay  on  the  list  a  long  time  (compared  to  the  fast  ones)  so  the 
collection  of  spikes  (usually  about  four)  representing  the  slow 
particles  is  not  changed  often.  In  Version  I,  the  instantaneous 
representation  of  the  low  velocity  peak  is  very  crude  (about  one 
spike  on  the  average)  but  the  entire  membership  of  the  list  is 
changed  after  eight  passages  of  the  test  particle.  Thus,  after 
four  complete  changes  of  the  list  (32  passages)  the  low-velocity 
peak  has  been  equally  well  represented  in  the  two  versions. 

Thus  it  is  difficult  to  choose  between  the  models,  and  the  choice 
must  depend  on  the  details  of  the  computation.  It  is  planned  to 
try  both  versions.  Version  II  may  possibly  be  preferable  because 
of  its  better  instantaneous  representation  of  the  low  velocity 
part  of  the  distribution. 
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